Nonlocal van der Waals functionals: The case of rare-gas dimers and solids 
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Recently, the nonlocal van der Waals (vdW) density functionals [M. Dion, H. Rydberg, E. 
Schroder, D. C. Langreth, and B. I. Lundqvist, Phys. Rev. Lett. 92, 246401 (2004)] have at- 
tracted considerable attention due to their good performance for systems where weak interactions 
are important. Since the physics of dispersion is included in these functionals, they are usually more 
accurate and show less erratic behavior than the semilocal and hybrid methods. In this work, several 
variants of the vdW functionals have been tested on rare-gas dimers (from He2 to KT2) and solids 
(Ne, Ar, and Kr) and their accuracy compared to standard semilocal approximations supplemented 
or not by an atom-pairwise dispersion correction [S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, 
J. Chem. Phys. 132, 154104 (2010)]. An analysis of the results in terms of energy decomposition 
is also provided. 
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I. INTRODUCTION 

Thanks to its relatively low cost/accuracy ratio, the 
Kohn-Sham (KS)i version of density functional theory 
(DFT) 2 - is the most used quantum method for the cal- 
culation of the geometrical and electronic properties of 
molecules, surfaces, and solids. The accuracy of the re- 
sults of a KS-DFT calculation depends primarily on the 
chosen approximation for the exchange-correlation func- 
tional (xc) E^c (see Ref. d for a recent review). Nowa- 
days, the most popular types of approximations for E KC 
are the semilocal [in particular the generalized gradient 
approximation (GGA)^] and hybrid functionals, 6 which 
very often give satisfactory results. However, it is well 
known that by construction none of these two approxi- 
mations account properly for the dispersion interactions, 
which arise due to the attraction between non-permanent 
dipoles, and that the results obtained with semilocal and 
hybrid methods on systems where dispersion interactions 
play a major role are often unreliable (see, e.g., Refs. 
andi). 

Therefore, efforts have been made to propose meth- 
ods within the framework of KS-DFT which explicitly 
account for the dispersion interactions (see Refs. 1914121 
for reviews). Among these methods, the simplest consist 
of adding to the KS-DFT total energy a dispersion term 
of the form 



disp 



= - E E ft mp (RA B ) 



a 



AB 



A<B n=6,8,10,. 



n-AB 



(1) 



where C^ B are the dispersion coefficients for the atom 
pair A and B separated by the distance Rab and f~ amp 
is a damping function preventing Eq. ([1]) to become too 
large at small Rab- The coefficients C^ B can be either 
precomputed (see, e.g., Refs. Iili - fl5l) or calculated using 
properties (e.g., electron density) of the system under 
consideration like in the exchange-hole dipole moment 
model (XDM) of Becke and Johnson^ 6 - or the method of 
Tkatchenko and Schemer.- 17 The DFT-D2 15 and DFT- 
D3i£ versions of Grimme are currently the most widely 
used of these methods. One of the advantages of most 



methods using Eq. ([TJ is to add a relatively negligible 
computational cost compared to the calculation of the 
KS-DFT energy. 

Another group of methods accounting explicitly of dis- 
persion interactions consist of adding a nonlocal term of 
the form 



E? = ~ 



p{r)<S> (r,r') p(r')d 3 rd 3 r' 



(2) 



to a LDA (local density approximation) or GGA corre- 
lation functional. In Eq. j2j, the kernel $ depends on 
quantities at r and r': 

4> (r,r') = 4> (p(r), p(r>), |Vp(r)| , |Vp(r')| , |r - r'|) . (3) 

The first functional of the form given by Eq. @, which 
could be applied to any type of systems was proposed 
by Dion et al. (DRSLL) ^ The DRSLL term was de- 
rived starting from the adiabatic connection-fluctuation 
dissipation theorem<22r— Originally, it was used in com- 
bination with revPBE 2 ^ (a reparametrization of the 
Perdew-Burke-Ernzerhof functional PBE 5 ) for exchange 
and LDA for correlation and the functional is named as 
vdW-DF in the literature (in Table fl] the composition 
of the functionals tested in the present work are given). 
vdW-DF was shown to be a clear improvement over the 
commonly used functionals, however it became also ob- 
vious that a serious shortcoming of vdW-DF is to sys- 
tematically overestimate the equilibrium distances . 24 ' 25 

Therefore, several attempts have been made to rem- 
edy this problem by combining the DRSLL nonlocal 
term with a more compatible semilocal functional or by 
proposing a new nonlocal term. For instance, Lee et 
al. (LMKLL) 26 proposed to modify slightly the DRSLL 
term (by changing the value of one parameter) and to 
use it in combination with PW86R 27 (a refitted version 
of the Perdew-Wang functional PW862&) . Their result- 
ing functional (called vdW-DF2) was shown to improve 
over the original vdW-DF. In Ref. H^, a new GGA ex- 
change functional (C09 x ) was proposed to be used with 
the DRSLL nonlocal term. The functional, C09 x -vdW, 
leads to more accurate results than vdW-DF for various 
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TABLE I. Composition of the tested exchange-correlation 
functionals. The VWN5 parametrization22 is used for LDA 
correlation. 



Functional 


Reference Exchange 


Correlation 


LDA 


1 and 33 


LDA 


LDA 


PBE 


o 


PBE 


PBE 


vdW-DF 


19 


revPBE 


LDA+DRSLL 


vdW-DF2 


26 


PW86R 


LDA+LMKLL 


C09 x -vdW 


29 


C09 x 


LDA+DRSLL 


optB88-vdW 


30 


optB88 


LDA+DRSLL 


RPBEc2/3+nl 


31 


RPBE 


iLDA+fPBE+DRSLL 


rVVIO 


34 


PW86R 


PBE+rVVIO 


PBE-D3 


18 


PBE 


PBE+D3 


revPBE-D3 


18 


revPBE 


PBE+D3 


B97D-D3 


18 


B97D 


B97D+D3 



types of systems. In Refs. |30| and l25l. Klimes el al. com- 
bined many GGA exchange functionals (already existing 
or newly proposed) with the DRSLL term. Particularly 
interesting are the functionals optB88-vdW and optPBE- 
vdW which lead to accurate results for both finite^ and 
extended^ systems. We also mention Ref. HH, where the 
nonlocal DRSLL term is used with the GGA RPBE ex- 
change functional 3 -^ and with either LDA for correlation 
or a linear combination of the LDA and PBE correla- 
tion functionals, the latter case leading to the functional 
named RPBEc2/3+nl (see Table U). 

Vydrov and Van Voorhis proposed their own nonlocal 
functionals, VV092 5 . and VVIO^ [also of the form given 
by Eq. @], which were constructed such that also the 
short-range regime of the van der Waals interactions is 
adequately described. VV10, when added to PW86R 27 
exchange and PBE correlation, has been shown to be par- 
ticularly accurate for finite systems (see, e.g., Ref. |37| ). 
However, its performance for solids is rather bad,— i 3 - 9 - 
therefore two parameters of the VV10 functional were 
modified such that the results for solids are improved. 39 

In this work, the results obtained for the equilibrium 
distance and interaction energy of rare-gas dimers and 
solids will be presented. The focus will be on the per- 
formance of several nonlocal functionals (listed in Table 
UJ with which only a few calculations on rare-gas sys- 
tems have been done up to now. The rare-gas systems 
are the prototypical van der Waals systems where the 
dispersion interactions are the only source of attraction 
between atoms and for which highly accurate ab initio or 
empirical results are available. The rare-gas dimers have 
been used numerous times for the testing of functionals 
for weak interactions (see, e.g., Refs. I4CH53I for extensive 
tests) , while tests on rare-gas solids are less common and 
recent^H 



II. METHODS 

The calculations were done with the Quickstep 
module 63 of the CP2K program package,— which is based 
on a mixed Gaussian and plane waves formalism.— More 
specifically, we used the Gaussian and augmented-plane- 
wave method (GAPW), 66 which allows for all-electron 
calculations. The calculations on the rare-gas dimers 
He2, Ne2, Ar 2 , and Kr 2 were done with the augmented 
correlation consistent polarized quintuple zeta (aug-cc- 
pV5Z) basis sets ; 67 ' 68 which lead to results very close to 
the basis set limit (see, e.g., Ref. f53h . In order to avoid 
the problem of linear dependence due to diffuse functions 
usually experienced in solids (as in the present case) , the 
calculations on solid Ne, Ar, and Kr were done without 
the augmentation functions by using the cc-pV5Z basis 
sets . 67 ' 68 The face-centered cubic (fee) structure was con- 
sidered for the rare-gas solids and we checked that using 
a unit cell comprising 32 atoms (2 x 2 x 2 of the fee four- 
atom unit cell) gives results which are very well converged 
with respect to the size of the supercell. 

The nonlocal term [Eq. @] was implemented accord- 
ing to the scheme of Roman-Perez and Soler^ 6 - 9 - which 
uses fast Fourier transforms, and therefore leads to cal- 
culations scaling as O(NlogN) (N is the number of 
points on the grid) instead of O (N 2 ^j for a direct evalu- 
ation of Eq. ([2]) in real space. Note that the method of 
Roman-Perez and Soler also leads to an efficient calcu- 
lation of the contribution of the nonlocal term to the 
KS-DFT potential (needed for the forces) and stress 
tensor* 6 " 9 -! 70 In our implementation, the nonlocal term 
is evaluated using only the smooth part of the electron 
density of the GAPW method. However, in Ref. H^, it 
was shown that within the projected-augmented wave 71 
method, plugin the all-electron density or the valence 
density into Eq. ([2]) leads to very similar results. Actu- 
ally, we checked that our results agree very closely with 
the results obtained with other codes when available [Ar 2 
with vdW-DF 1 - 9 - 3 ^^ and (r)VV10 34 < 36 and Kr 2 with 
vdW-DF-i 9 - 7 ^ and W1<M]. 

In addition to the functionals already introduced in 
Sec. HI we also mention rVVlO^ which is a revised 
version of VVIQ 3 ^, such that its evaluation can also be 
done with the method of Roman-Perez and Soler. It was 
shown (Ref. l34| ) that VV10 and rVVIO give very sim- 
ilar results. rVVIO is among the functionals tested in 
the present work (Table [I]). For comparison purposes, we 
also considered the standard functionals LDA and PBE, 5 
as well as the dispersion-corrected functionals PBE-D3, 
revPBE-D3, and B97D-D3 (B97D 15 i s a reparametriza- 
tion of B97— ), where D3 refers to the third set of param- 
eters C£ B [in Eq. (JlJ] proposed by Grimme (the three- 
body term was included in our calculations) .— Note that 
we used the VWN5 parametrization 3 - 3 - for the LDA cor- 
relation. Finally, we mention that LIBXC, a library of 
exchange-correlation functionals, 7 " 6 - has been used for the 
evaluation of some of the semilocal functionals in Table 
III 
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III. RESULTS 
A. Rare-gas dimers 

The interaction energy curves of the rare-gas dimers 
He2, Ne2, Ar 2 , and Kr 2 are displayed in Fig. Q]and the 
corresponding values at the minimum (equilibrium dis- 
tance Ro and binding energy AE) are shown in Table 
HJ The KS-DFT results are compared to very accurate 
reference (theoretical or experimental, see Ref. [z3 for 
details) results. 

Discussing first the results obtained with the 
(semi)local functionals, it is already knowni^ that LDA 
strongly underestimates the bond lengths and overesti- 
mates the binding energies of all dimers, this trend be- 
ing systematically observed with LDA for intermolecular 
complexes. For He2 and Ne2, LDA leads to binding en- 
ergies which are one order of magnitude too large and 
to distances which are about 0.5 A too small. Among 
the countless semilocal and hybrid functionals tested on 
rare-gas dimers, PBE (a GGA free of any empirical pa- 
rameter) is one of the most accurate (or least inaccu- 
rate, see Refs. HUES andl47l for extensive tests). Still, 
PBE accuracy can not be considered as satisfying, Ne 2 
excepted, since it largely overbinds He2 and underbinds 
Ar 2 and Kr2. In the group of semilocal and hybrid func- 
tionals (which do not include the physics of dispersion in- 
teractions), it was showr4£ that the hybrid B97-1— and 
meta-GGA hybrid M05-2X 79 arc also among the best for 
rare-gas dimers, but as in the case of PBE the results are 
in some cases rather inaccurate. Finally, from all pre- 
vious studies on rare-gas dimers, we can conclude that 
there is no semilocal or hybrid functional that can be 
considered as reliable. 

Turning now to the results obtained with the six non- 
local functionals, we can see that a large range of re- 
sults can be obtained. vdW-DF largely overbinds all four 
dimers, and while the bond length is reasonable for HC2 
and Ne2, it is too large by 0.2—0.3 A for Ar2 and Kr2. 
vdW-DF2 improves over vdW-DF for AE by reducing 
the overbinding by a factor of two, but now the bond 
lengths of He 2 and Ne2 are clearly too short. The bond 
lengths obtained with the C09 x -vdW functional are as 
inaccurate as the LDA ones, but with the opposite trend 
(overestimation ranging from 0.2 A for He 2 to 0.7 A for 
K1-2). The C09 x -vdW interaction energies are not par- 
ticularly accurate except for Ar 2 . optB88-vdW leads to 
quite accurate results for the binding energy AE, how- 
ever, the bond lengths Rq are too large (in particular for 
He 2 with 0.5 A of error). The binding energies obtained 
with the RPBEc2/3+nl functional constitute a disaster 
since they are even larger than LDA values, while the 
equilibrium bond length is accurate for Ar2 and Kr2, but 
not for the two lighter dimers. Among all tested func- 
tionals in this work, rVVIO is clearly the most accurate 
one. From Fig. [lja), we can see that for He2 the rVVIO 
and reference curves coincide very closely along the whole 
range of considered intermolecular distances and corre- 



spond to the same binding energy (0.9 meV). For the 
other dimers, also both the bond lengths and interaction 
energies are very accurate. The largest error in AE is 
only 2 meV (for Ne2 and KX2). 

Concerning the three DFT-D3 methods that we consid- 
ered, revPBE-D3 leads to a rather accurate bond length 
for He2, but overestimates AE. For the other three 
dimers, revPBE-D3 yields values for Rq, which are too 
large by 0.1—0.2 A, but quite accurate values for AE. 
Overall, PBE-D3 leads to values which are less satisfying 
than revPBE-D3. B97D-D3 leads to results which are 
quite similar to revPBE-D3, but overestimates the bond 
lengths even more for Ne2 and Ar2. When compared 
to the nonlocal functionals, revPBE-D3 and B97D-D3 
seem to show more stability in the results, except when 
compared to rVVIO which is by far the most accurate 
functional. 

Among the previously published works on the testing 
of DFT functionals on rare-gas dimers, we should men- 
tion the results obtained by Kannemann and Becke 50 
with the functional PW86xPBEc-BJ (results also shown 
in Table HI| . where BJ refers to the XDM model for dis- 
persion of Becke and JohnsonJ£ Their calculated bond 
lengths and binding energies are in very close agreement 
with the reference results, and actually the accuracy of 
rVVIO and PW86xPBEc-B J can be considered as similar. 
However, it is important to note that the two adjustable 
parameters in the PW86xPBEc-BJ functional were de- 
termined by minimizing the error for AE of a set of ten 
rare-gas dimers (all combinations involving He, Ne, Ar, 
and Kr). rVVIO (VV10) also contains two parameters, 
but one of them was adjusted such that the mean error 
of Cq A coefficients for a set of 54 species (among them 
He, Ne, Ar, and Kr) is minimized,^ while the other was 
determined using the S22 set of noncovalent complexes^ 
which does not contain any rare-gas atoms. Therefore, 
(r)VV10 was certainly not adjusted exclusively on rare- 
gas systems, which makes its excellent performances on 
these systems even more impressive. 



B. Rare-gas solids 

The results for the rare-gas solids Ne, Ar, and Kr are 
shown in Fig. [5] and Table Hill The reference results were 
obtained from coupled cluster with single, double, and 
perturbative triple excitations [CCSD(T)] calculations. 81 
For a meaningful comparison of our KS-DFT results with 
the CCSD(T) results, the zero-point energy calculated in 
Ref. l8l| has been removed from the CCSD(T) results. 

As for the dimers, LDA leads to severe underestima- 
tion of the lattice constant ao and overestimation of the 
cohesive energy AE for all three solids. Note, however, 
that in some cases LDA can, at a qualitative level, give 
relatively correct results for solids which are bound by 
weak interactions. Such examples include layered solids 
like graphite (see, e.g., Refs. |H and [55|). As already 
observed in Ref. [5f| PBE gives essentially the same co- 
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TABLE II. Equilibrium bond length Rq (in A) and interaction energy AE (in meV and with opposite sign) of rare-gas dimers 
calculated from various functionals and compared to accurate reference values and results obtained from the exchange-hole 
dipole moment model of Becke and Johnson (BJ). 



He 2 Ne 2 Ar 2 Kr 2 



Functional 


Ro 


AE 


Ro 


AE 


-Ro 


AE 


Ro 


AE 


LDA 


2.40 


9.6 


2.64 


20.4 


3.40 


30.9 


3.68 


36.7 


PBE 


2.76 


3.2 


3.08 


5.6 


4.00 


6.3 


4.36 


6.9 


vdW-DF 


2.82 


6.6 


3.07 


14.1 


3.92 


23.1 


4.27 


26.2 


vdW-DF2 


2.75 


2.8 


2.95 


9.2 


3.75 


18.3 


4.09 


22.3 


C09 x -vdW 


3.19 


4.1 


3.51 


6.6 


4.37 


11.5 


4.71 


13.4 


optB88-vdW 


3.48 


0.5 


3.30 


3.0 


3.93 


11.7 


4.20 


16.1 


RPBEc2/3+nl 


2.66 


11.2 


2.93 


23.1 


3.77 


34.4 


4.10 


38.0 


rVVIO 


2.92 


0.9 


3.01 


5.6 


3.73 


13.9 


4.00 


19.3 


PBE-D3 


2.66 


5.7 


3.01 


9.9 


3.88 


15.3 


4.16 


19.3 


revPBE-D3 


2.90 


3.0 


3.20 


5.6 


3.93 


12.8 


4.18 


17.9 


B97D-D3 


3.01 


2.4 


3.33 


4.3 


3.99 


11.3 


4.18 


17.2 


PW86xPBEc-BP 


3.01 


0.8 


3.12 


3.8 


3.84 


11.2 


4.07 


17.0 


Reference^ 


2.97 


0.9 


3.09 


3.6 


3.76 


12.4 


4.01 


17.4 


a Reference |5& 
b Reference [73. 



TABLE III. Equilibrium lattice constant ao (in A) and cohesive energy AE (in meV/atom and with opposite sign) of rare- 
gas solids calculated from various functionals and compared to reference [CCSD(T)] results as well as the PBE-TS and RPA 
methods. 



Ne Ar Kr 



Functional 


ao 


AE 


ao 


AE 


ao 


AE 


LDA 


3.86 


92 


4.94 


136 


5.36 


164 


PBE 


4.55 


25 


5.93 


25 


6.42 


25 


vdW-DF 


4.32 


101 


5.49 


163 


5.96 


184 


vdW-DF2 


4.17 


65 


5.29 


130 


5.75 


157 


C09 x -vdW 


4.90 


51 


6.00 


83 


6.39 


96 


optB88-vdW 


4.24 


59 


5.24 


143 


5.63 


181 


RPBEc2/3+nl 


4.19 


146 


5.35 


222 


5.80 


246 


rVVIO 


4.19 


49 


5.17 


117 


5.53 


162 


PBE-D3 


4.37 


53 


5.58 


84 


5.93 


108 


revPBE-D3 


4.66 


32 


5.62 


71 


5.89 


104 


B97D-D3 


4.78 


26 


5.69 


66 


5.87 


104 


PBE-TS^ 


4.42 


43 


5.51 


83 


5.90 


97 


RPA(PBE)^ 


4.5 


17 


5.3 


83 


5.7 


112 


CCSD(T£ 


4.297 


26 


5.251 


88 


5.598 


122 



a Reference 1 6 ll. 

b RPA ener gy evaluated with PBE orbitals and eigenvalues— 
c Reference 1 8 ll. 



hesive energy (25 meV/atom) for the three solids (very 
large underestimation for Ar and Kr). This is somewhat 
similar to what is observed for the corresponding dimers 
(see Table HI)) . The PBE lattice constants are by far too 
large (by more than 0.7 A for Ar and Kr). 

Concerning the nonlocal functionals based on the 
DRSLL or LMKLL kernels, the observations are the fol- 



lowing. vdW-DF and vdW-DF2 clearly overbind the 
rare-gas solids, C09 x -vdW totally fails for the lattice con- 
stant, and RPBEc2/3+nl leads to the largest overbinding 
(as for the dimers). optB88-vdW leads to quite accu- 
rate values for ao , but overestimates the cohesive energy 
rather strongly (50%— 100%), while in the case of the 
dimers optB88-vdW was quite good for the interaction 
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FIG. 1. Interaction energy curves for (a) He2, (b) Ne2, (c) Ar2, and (d) Kr2 obtained from various functionals and compared 
to reference results^ (black line without symbols). 



energy. 

The rVVIO nonlocal functional seems to be again supe- 
rior to the other nonlocal functionals. The bond lengths 
are rather good (albeit too short by ~ 0.1 A for Ne), while 
the cohesive energies are too large for all three solids, but 
the error is smaller than for the other nonlocal functionals 
except C09 x -vdW. The cohesive energies obtained with 
the DFT-D3 methods are quite accurate (PBE-D3 for Nc 
excepted), however the lattice constants are consistently 
too large by more than 0.3 A in most cases. 

Also shown in Table |TT] are the results from Ref. [6l| 
obtained with the Tkatchenko and SchefHer— (TS) ap- 
proach using PBE for the semilocal part. We can see 
that the PBE-TS results are similar to the results from 



PBE-D3 for both the lattice constant and cohesive en- 
ergy. For completeness, we also show in Table [TT| the 
values from the non-DFT method RPA (random-phase 
approximation).™ The RPA bond lengths are somehow 
overestimated, but the cohesive energies are very close 
to the CCSD(T) values. The RPA method is superior 
to the KS-DFT methods considered in the present work, 
but leads to calculations which are obviously much more 
expensive. Finally, we mention the DFT+XDM results 
from Ref. [H for the lattice constant ao, where the XDM 
dispersion correction was added to two different GGA 
functionals. The results are good only for Kr, whereas 
for Ne and Ar rather inaccurate values were obtained. 
As a summary of the results on rare-gas solids, the 
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FIG. 2. Cohesive energy curves for (a) Ne, (b) Ar, and (c) Kr obtained from various functionals and compared to CCSD(T) 
results*- (black square). 



rVVIO nonlocal functional seems to be a relatively good 
choice (at least compared to the other functionals), but 
leads to non-negligible overestimations of the cohesive 
energy, and this more than in the case of the rare-gas 
dimers. 



IV. FURTHER DISCUSSION 
A. Energy decomposition 

In order to have more insight into the results of Sees. 
IIII Al and IIII B[ we now consider the various contri- 
butions to the interaction energy AE and their rel- 
ative importance. Figure [3] shows for the Ar dimer 
(similar trends are observed for the other rare-gas 
dimers and solids) the contributions to AE coming from 
the (semi)local exchange-correlation functional (AE^ C ), 
the (atom-pairwise or nonlocal) dispersion energy term 
(A££ is P), the sum of these two (AE XC = AE^ c +AEf sp ), 
the rest of the terms (AE Icst ) representing the kinetic 
and electrostatic energies, and the sum AE Icst + AE^ C . 
The sum of AE XC and AE Ies t gives the total interaction 
energy AE shown in Fig. HJc). 

Around the experimental equilibrium bond length (~ 
3.8 A) all terms seem to be of roughly equal importance 
and the absolute values range from 10 to 50 meV de- 
pending on the functional. However, for smaller bond 
lengths R, the curves AE^ C and AE^t vary faster than 
A£^ lsp and these terms become much more important. 
For instance, at R — 2.0 A (not shown), the magnitude 
of AE ICBt and AE*[ is around 8000 and 4000 meV, re- 
spectively, while for AE^ lsp it is smaller than 5 meV for 



the DFT-D3 methods and between 100 and 300 meV for 
the nonlocal methods. Note that the different behavior 
of AE^ is P for the DFT-D3 method at small values of R 
is due to the damping function /^ amp in Eq. (p}. 

In Fig. [2Jb), we can see that among the nonlocal 
functionals, rVVIO and DRSLL lead to the smallest and 
largest values (in magnitude) for AE^ lsp , respectively. In 
the range of intermolecular distances that we considered, 
AE^ isp is negative, however, the energy (the component 
of the total energy) is always positive in the case of the 
nonlocal functionals [Eq. J5J], while the values are nega- 
tive for the atom-pairwise DFT-D3 method [Eq. JTJ]. 

The sum of the terms AE^ C and AE rest , which rep- 
resents the interaction energy of Ar2 calculated without 
the dispersion term, is shown in Fig. [He). Among the 
scmilocal functionals only PBE and RPBEc2/3 yield rea- 
sonable interaction energies, while all other functionals, 
except LDA barely bind or do not bind at all the two 
Ar atoms. Actually, it is clear that in order to avoid 
overbinding (due to double counting) when using a dis- 
persion term in the total-energy expression, it should 
be combined with a semilocal functional which leads to 
strongly underestimated interaction energyJ^ 

From Fig. |31 we can also infer that the differ- 
ences in AE between the functionals can not be un- 
derstood by looking exclusively at the contribution from 
the exchange-correlation energy. Indeed, the curves for 
Ai5 rest (whose analytical form is the same for all func- 
tionals) show differences which are as strong as for AE^ C 
and AE^ lsp . Actually, the differences in the Ai? rost 
curves are a reflection of the corresponding exchange- 
correlation potentials u xc = SE XC /Sp used in the KS- 
DFT equations. It is known that for the calculation of 
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FIG. 3. Contributions to the Ar2 interaction energy coming 
from (a) the semilocal (si) xc functional, (b) the dispersion 
energy [Eqs. (JXJ) or @], (c) the total (si plus dispersion) xc 
functional, (d) the rest (kinetic plus electrostatic), and (e) the 
sum of si and rest. The addition of (c) and (d) gives the total 
interaction energy of Fig. [T^c) . 



FIG. 4. Same as Fig. but obtained from non-self-consistent 
calculations by plugin the PBE orbitals and electron density 
into the functionals. 



properties depending on total energies, the results usu- 
ally do not depend too sensitively on the orbitals and 
electron density plugged into the total-energy functional 
(but one needs to be very careful with this statement). 
However, the individual components of the total energy 
show much stronger sensitivity, but these variations tend 
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FIG. 5. Three-body interaction energy of Ar3 (equilateral 
triangular configuration) plotted against the Ar-Ar distance. 
The CCSD(T) results are from Ref. [H. ATM (black line 
without symbols) is the Axilrod-Teller-Muto term given by 
Eq. ©. ' 



to cancel among the different terms. For A1-2, we also 
performed non-self-consistent calculations by evaluating 
all functionals with the PBE orbitals and electron den- 
sity (results shown in Fig. 2]). The resulting equilib- 
rium bond lengths and binding energies are essentially 
the same as their self-consistent counterparts. However, 
from Fig. Eta) we can see that for some functionals 
(LDA and C09 x -vdW in particular), the AE^ C curve 
is quite different to the one obtained self-consistently 
[Fig. Eta)]. For the nonlocal dispersion terms (DRSLL, 
LMKLL, and rVVIO), basically no difference between the 
self-consistent and non-self-consistent calculations can be 
seen, which is maybe due to the fact that these terms are 
evaluated only with the smooth part of the electron den- 
sity. 



B. Three-body interaction energy 

The leading term in the many-body contribution to the 
interaction energy is the three-body nonadditive energy 
AE3 . In the case of a simple trimer this term is calculated 
as the atomization energy of the trimer minus the sum of 
the atomization energies of the three dimers. If the three 
atoms in the trimer are identical, then 



The asymptotic behavior of the dispersion component 
of AE 3 is given by the Axilrod-Teller-Muto^s (ATM) 
triple-dipole term 



AE. 



ATM 



_ r<ABC 
— °9 



1 



COS 



(Oabc) cos{9 B ca) cos(6 C ab) 



{RabRbcRcaY 



AE* = E, 



trimer 



dimcr 



(4) 



(5) 

where Oijk and Rij are the angles and side lengths of the 
triangle formed by the trimer and C$ BC is the triple- 
dipole constant. 

Using Eq. (j4}, we calculated AE3 for the Ar trimer 
at equilateral geometry, and in Fig. [5] the results are 
compared to the accurate CCSD(T) values from Ref. 
l83l as well as the asymptotic ATM term [Eq. ([5])] with 
Cq BC — 521.7 au4& In general, the three-body interac- 
tion energy AE3 of a trimer in this configuration is the 
largest contribution to the many-body cohesive energy 
of the corresponding solid in the fee structure. For in- 
teratomic distances R larger than ~ 3.5 A, we can see 
that most functionals strongly overestimate (too positive 
values) the three-body energy AE$. The exceptions are 
optB88-vdW which leads to negative values for the whole 
range of interatomic distances R that we considered and 
rVVIO which seems to be the best of the considered func- 
tionals. Closely around the equilibrium interatomic dis- 
tances in the dimer and solid (~ 3.7-3.75 A), the rVVIO 
values are close to the CCSD(T) and ATM values. How- 
ever, the maximum of the AE3 curve is at ~ 3.2 A for 
rVVIO (and overestimated) , while it is at ~ 3.7 A for 
CCSD(T). In Refs. l57l andl59i. other functionals were con- 
sidered for the calculation of AE 3 , but none of them lead 
to results in qualitative agreement with the CCSD(T) re- 
sults. 



V. SUMMARY 

We have presented the results of KS-DFT calculations 
on rare-gas dimers and solids. The focus was on the 
performance of nonlocal vdW functionals for the equi- 
librium bond length and binding energy. The main con- 
clusions are (a) overall the rVVIO functional is the one 
performing the best, (b) some others (e.g., C09 x -vdW 
or RPBEc2/3+nl) can perform very badly, and (c) the 
considered DFT-D3 methods show good accuracy for the 
interaction energy, but seem to lead to some (slight) over- 
estimation of the bond lengths. In Sec. IIV A[ we pre- 
sented an analysis by decomposing the interaction energy 
into its components in order to estimate the relative im- 
portance of each term, and in Sec. IIVBI it was shown 
that rVVIO gives also reasonable values (at least close 
to the equilibrium geometry) for the three-body nonad- 
ditive energy, while the other DFT functionals are very 
inaccurate. 

Considering the results on rare-gas systems obtained in 
the present work and from previously published papers, 
the (r)VV10 nonlocal functional seems to be the most ac- 
curate among the DFT methods. It was already shown to 
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be accurate for many other finite systems^ while for lay- 
ered solids it is necessary to modify its parameters (and 
eventually combine it with another semilocal functional) 
to get accurate results^ 
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